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Summary 

An efficient implicit finite-difference algorithm for the gasdynamic equa- 
tions utilizing matrix reduction techniques is presented. A significant reduc- 
tion in arithmetic operations is achieved while maintaining the same favorable 
stability characteristics and generality found in the Beam and Warming ap- 
proximate factorization algorithm. Steady-state solutions to the conservative 
Euler equations in generalized coordinates are obtained for transonic flows 
about a NACA 0012 airfoil. The theoretical extension of the matrix reduc- 
tion technique to the full Navier-Stokes equations in Cartesian coordinates 
is presented in detail. Linear stability, using a Fourier stability analysis, is 
demonstrated and discussed for the one-dimensional Euler equations. It is 
shown that the method offers advantages over the conventional Beam and 
Warming scheme and can retrofit existing Beam and Warming codes with 
minimal effort. 

I. Introduction 

Numerical algorithms for solving the compressible Euler and Navier-Stokes 
equations have received considerable attention over the last decade. Both 
implicit and explicit time advance and iterative techniques have been utilized. 
Explicit methods typically have lower arithmetic operation counts but more 
restrictive stability bounds than implicit methods. Use of implicit methods 
is generally advantageous for flow situations in which an explicit time step 
limitation would be much less than the time step necessary for the desired 
accuracy. This situation frequently arises in the solution of unsteady flows 
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that are characterized by low frequency motion, boundary condition forcing, 
and high Reynolds number viscous effects. Implicit algorithms can often be 
readily adapted to be efficient relaxation schemes for steady state applications 
as well. 

The purpose of this paper js to present a technique that substantially 
reduces the arithmetic operation count, but otherwise retains the stability 
characteristics of a Beam- Warming approximate factorization implicit al- 
gorithm for the Euler equations. This new method extends a matrix splitting 
of Steger (ref. 1), previously restricted to Cartesian coordinates, such that 
equation decoupling is achieved for the conservative form of the Euler equa- 
tions. The equation decoupling (or matrix reducibility) results in a significant 
improvement in the efficiency of the algorithm. By utilizing a local trans- 
formation, the equation decoupling techniques are extended to the Euler 
equations in generalized coordinates. Computation of transonic flow about a 
NACA 0012 airfoil is used to verify that no numerical stability is lost with the 
new reducible implicit algorithm. Extension of the technique to the Navier- 
Stokes equations is also indicated in Appendix A. 

This work was partially funded by the Air Force Office of Scientific 
Research under AFOSR-82-0254. The discussions and support of Drs. Joseph 
Shang and Wilbur Hankey of the Air Force Flight Dynamics Laboratory are 
also gratefully acknowleged. 

II. Background 

In this section we review a sound speed and velocity splitting concept 
for an approximate factorization implicit algorithm in Cartesian coordinates 
and describe how it can be used to improve computational efficiency. The 
extension for general coordinate systems is described in the following section. 


A. Beam and Warming Algorithm 

The conservative form of the Euler equations for a perfect gas can be 
written in Cartesian coordinates as 

d t Q + d x E + d y F = 0 (2.1) 


where 
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and 


P = (l— l)(e - | p(u 2 + z/ 2 )j (2.3) 

Here p is the fluid density, u and v are the Cartesian velocity components, e 
is the total energy per volume, and p is the pressure. 

A general-purpose implicit finite difference scheme (ref. 2 and 3) for equa- 
tion (2.1) that can be readily adapted to either steady or unsteady flow ap- 
plications is given by 


I + h6 x A n + 

-At 


I + h8 y B n + D& 


A Q n 


6 x E n + 8 y F n + D^Q n 


(2.4) 


where 

£>< 4 > = e e At[(V Ag + (VA$ (2.5a) 

and 

D® = -6iAt(VA) x , D® = -eiAt{VA) y (2.5b) 

Here 8 X and 5 y are three-point, second-order accurate, central space difference 
operators, while A and V denote forward and backward space differences. 
Either first- or second-order time accurate differencing may be used with 
h = At or At/2; alternately, the time step may be used as a relaxation 
parameter and may vary in space as well. With the use of central space 
differencing numerical dissipation is needed and implemented with the second 
and fourth order differences given above. The parameters e e and control 
the amount of numerical dissipation. These parameters may be constants or 
may vary with the solution gradients (in which case they are moved inside 
the operators so as to maintain conservative form). Steady state convergence 
can be greatly enhanced by more implicit treatment of the dissipation terms 
and by use of space- varying scaling parameters. 

In equation (2.4), local time linearizations have been utilized to avoid 
solving nonlinear equations between time or iterative levels n to n-j-1 by 
expanding E and F in terms of A Q as 


E n+ 1 = E n + A n {Q n+1 - Q n ), A = dE/dQ (2.6a) 

F n+1 = F n + B n (Q n+1 - Q n ), B = dF/dQ (2.6b) 

The left hand side operators of equation (2.4) form a linear system of 
equations that must be solved at each time or iteration level. Although 
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the matrix band structure of this system has been altered by approximately 
factoring into one-dimensional-like operators, the inversion process has been 
simplified. The approximate factorization (AF) reduces the inversion work 
but sometimes at the cost of limiting the time step At because of either 
time accuracy considerations or because of reduced inefficiency in relaxation 
applications. 

In practice the difference equations (2.4) are solved in the ADI algorithm 
form 


L X AQ* = RHS{ 2.4) (2.7a) 

L y AQ n = AQ* (2.7b) 

Variants of this algorithm include higher-order difference approximations 
(including pseudo-spectral right-hand-side operators (ref. 4)), space varying 
At for steady state applications, addition of viscous terms, etc. [ref. 5-9] 

B. Diagonalization 

Efforts have been made to reduce the inversion work over and above 
what is obtained with approximate factorization. In particular, Pulliam 
and Chaussee (ref. 7) have used eigenvector-similarity transforms to reduce 
the block tridiagonal matrices to scalar tridiagonals (see also ref. 8-10). 
Specifically, the eigenvectors of the A and B Jacobian matrices, denoted as X 
and Y, are further approximately factored into the left-hand- side of equation 
(2.4) as 
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8 x E n + 8 y F n -\-D^Q n 
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where 


A a =X~ 1 AX and A b = Y~ 1 BY 


— i 


As fig. 1 illustrates, the number of operations in solving a block tridiagonal 
(that uses only L-U decomposition) increases with the cube of the block size, 
a formula given in ref. 11 indicates the work increases as -^m 3 -j- |m 2 — £m, 
where m is the dimension of the block. Specifically our own coded 4 x 4 
block solver requires 370 operations per entry (i.e., grid point) while the lxl 
tridiagonal needed in equation (2.8) requires only 9 operations per entry (or 38 
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operations for 4 scalar tridiagonals). Although the eigenvector matrices must 
be formed and multiplied through (each 4x4 multiply requires 28 operations 
per entry), the arithmetic operation count of equation (2.8) is considerably 
reduced from that of equation (2.4). For block tridiagonal matrices subject 
to periodicity conditions and for block pentadiagonal matrices, the saving is 
even more significant. 

C, Sound Speed and Velocity Splitting 

An alternate method for reducing the inversion (i.e., solution) work was 
proposed and demonstrated in ref. 1. This method was also motivated by 
similarity transforms, but it does not explicitly use them in the algorithm. 
Instead, it involves splitting the Jacobian matrices A and B into velocity and 


sound speed (or pressure parts) as 

.A — A u “f- A c (2.9a) 

B = B v -j- B c (2.9b) 

such that the eigenvalues a(A) and o(B) are 

o{A) = a(A u ) + a(A c ) (2.10a) 

<j(B) = o(B.) + a(B c ) (2.10b) 

with 

a{A u ) — {u, u, u,u), a{A c ) = (0, c,0, — c) (2.11a) 

a{B v ) — (v,v, v, v), cr(B c ) = (0, 0, c, — c) (2.11b) 


The matrices A U) A C ,B V , and B c were found from a natural reduced form of 
the equations in nonconservative form. Specifically A c and B c were split, as 
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*41 = M« 2 + y2 )/2] - 1 v P/[p(l - l) 2 ] 

^43 = 7P/[p(7- l) 2 ]- V 2 

while A u and B v are 

A u — A— A c (2.13a) 

B V = B-B C (2.13b) 

This splitting produces matrices A u and B v that are more complex than 
A and B. However, Steger (ref. 1) noted that Q is an eigenvector of both 
matrices, i.e., 

A^Q == uQ 
BvQ — vQ 

which prompted the ad hoc substitution 

AQ — (ul -f- A C )Q (2.15a) 

BQ = (vI + B c )Q (2.15b) 

The matrices ul + A c and vl + B c have a reduced form that simplifies 
inversion compared to A and B. 

Insertion of equation (2.15) into the equations for local linearization of the 
Jacobians, equation (2.6), produces 


(2.14a) 

(2.14b) 


E n+1 = E n + ( u/ _j_ A c ) n {Q n+1 - Q n ) (2.16a) 

F n+ 1 = F n + (vl + B c ) n (Q n+1 - Q n ) (2.16b) 

Utilizing these linearizations in the basic, algorithm equation (2.4) 

L x L y AQ — RHS(2A) (2.17) 

gives new left- hand- side operators L x and L y that are easier to solve 
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(2.18b) 
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The end result of this splitting is that the new operators L x and L y form 
matrices that no longer require 4x4 block tridiagonal inversions. Leaving 
out the dissipation terms to illustrate the structure of these operators, we 

obtain 
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(2.19a) 


(2.19b) 


where a c and b c are the respective elements of A c and B c given by equation 

(2.12). 

For the L x operator, for example, the first and third rows decouple from the 
system and can be solved as scalar tridiagonal matrices with their respective 
right-hand-sides. Once these rows are solved, the elements of the first and 
third columns can be moved to the right-hand-side. The second and fourth 
equations remain coupled and are solved as a 2 x 2 block tridiagonal matrix. 
The dissipation terms also form diagonal tridiagonals and so do not alter the 
structure. The block decoupling of the L y operator is even more conspicuous 
and is inverted (i.e., solved for) in a similar manner. 

Substitution of the left-hand-side operators given by equation (2.18) in 
place of those of equation (2.7) results in a substantial reduction in arithmetic 
operations. Onr typical 2x2 block tridiagonal requires 55 operations per 
point, so the overall inversion, including the two scalar tridiagonals, requires 
73 operations per entry. Because the two scalar tridiagonals have identical 
coefficients, this work can be even further cut by solving them together. 

The matrix splitting produces the flux vectors 
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Note that we do not reproduce the matrix splitting in equation (2.9) (with 
A c and B c defined from equation (2.12)) by taking the Jacobians of equa- 
tions (2.21). Surprisingly, linear stability analysis (Appendix B) as well as 
numerical experiment have shown that the use of the Jacobian matrices for 
A c and B c is unsatisfactory. However, the same stability analysis and numeri- 
cal experimentation have confirmed that no numerical stability degradation 
occurs by using the substitute linearization matrices, equation (2.15). But 
the linearization in equation (2.16) is only first-order accurate because of the 
substitution of equation (2.14). 

The matrix splitting concept can be extended to three dimensions. In this 
case the difference equations can be represented by 


L x L y L z (Q n+1 - Q n ) = - A t[8 x E n + 6 y F n + 6 z G n + D^Q n ] (2.21) 

Associated with each operator are 5 X 5 block tridiagonal matrices, with each 
matrix inversion requiring about 700 operations per grid point. Typically the 
inversion represents 70 to 80 per cent of the overall work per point. Replacing 
each Jacobian A by ul -j- A C) etc., reduces the inversion cost from 700 to 
82 operations per entry of each block tridiagonal. By combining the scalar 
tridiagonals, this can be cut to 74 operations. 
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TFT . Sound Speed and Velocity Splitting in Steady Generalised Coordinates 

The two-dimensional Euler equations can be transformed from Cartesian 
coordinates to steady curvilinear coordinates using new independent variables 

T = t 

Z = Z(x,y) (3.1) 

n = 7j{x,y). 

The Euler equations written in steady generalized curvilinear coordinates 
are 

d T Q + df,E + djjF = 0 (3.2) 
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where J is the Jacobian £ x 7) y — £, y rj x and 

U — + V = rj x u + 7] y v 


are unsealed contravariant velocity components. 

Application of the AF implicit scheme, equation (2.4), to equation (3.2) is 
given by 
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and 
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The Jacobian matrices of E and F can be given in terms of the Jacobian 
matrices A and B as 

A = Z x A+Z y B (3.4a) 
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B=7] x A+ri y B (3.4b) 

Making the substitutions of equation (2.15 ) into AQ and BQ results in 

AQ = (i UI + ( X A C + ( y B c )Q (3.5a) 

BQ = (VI + r)xA c + Vy B c )Q (3.5b) 


Because of the linear combination of both A c and B c being present in 
equation (3.5), only 3x3 reducibility is achieved. When applied in the 
implicit AF algorithm, this results in a savings as indicated in ref.l, but not 
as great a saving as was achieved in Cartesian coordinates. Moreover, in 
three dimensions one can only reduce to a 4 X 4 block tridiagonal. 

The transformed equations given above have only been transformed in their 
independent variables. That is, Cartesian velocity or momentum components 
have been kept, and, as a result, the so-called strong conservation law form is 
maintained. Suppose, for example, we had used orthogonal body coordinates 
s and n, and that we had transformed the momentum components as well. 
We would then obtain equations with coordinate source terms, but otherwise 
the equations would look like their Cartesian counterparts and the Jacobian 
matrices could be reduced as before. We can achieve this same effect with 
equation (3.2) by multiplying through by the matrix 


T 0 0 0- 

o (x o 

0 J] x rjy 0 

-0 0 0 1 - 


(3.6) 


This matrix transforms u,v momentum components to U,V components. 
With multiplication of equation (3.2) by C we obtain 
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with 

H = -C\(C i r'E + (C^)~ l F] 

These equations look much like their Cartesian counterparts except for 
the coordinate source term and the appearance of two extra pressure terms. 
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However, for orthogonal coordinates, V£ • V?? = 0, and the extra pressure 
terms then drop out. Although we will not give the details here, one can 
anticipate from the form of the equations that the flux Jacobians in the 
implicit algorithm are 2x2 reducible (here the Jacobian matrices are formed 
with respect to pU and pV components, not pu and pv components). So if 
the coordinate generated source term is either treated explicitly or properly 
distributed in the or operators, the inversion efficiency of the previous 
section can be obtained. 


Generally the transformed coordinates will not be orthogonal and multi- 
plication through by C will not lead to the 2x2 reduced matrix operators. 
However, if we use a transform matrix that brings out a mixture of con- 
travariant and covariant velocity components, the extra pressure terms will 
also disappear. Although it may not be immediate^ obvious, we can again 
obtain the reduced matrix operator structure without requiring the coor- 
dinates to be orthogonal. 

Consider the transform matrix C given by 
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where 
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t 2 = y/vt) • Vt? 


and the inner 2 x 2 of C transforms u and v into the covariant velocity 

components 


u = \rjyU - Vxv), V — i x '{-ZyU + txV) 

Since the determinant of C is Jl 2 1 , its inverse exists for all mappings of 

interest and is 
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Multiplying equation (3.2) by C we obtain 

d T Q -f d<:E + d n F = II (3.10) 

where 

// == + C~ X F) 

and Q — CQ, E = CE, F — CF. The flux vectors written out are 
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Unlike equation (3.7), extra pressure terras are not picked up in the momen- 
tum equations. 

The previous numerical algorithm extended to equation (3.10) is given by 
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where we have lagged the source term and where A = dEjdQ and B — 
dF/dQ. 

It can be verified that A — CAC 1 and B — CBC 1 . Using these 
relations as well as equation (3.5) we find 


AQ = ( UI + A C )Q BQ = (VI + B C )Q 


(3.13) 
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where 


“41 = u\tf + “ 2 )/2 - ipMi - I) 2 )] 

“42 = khpClMl ~ l) 2 ) - U 2 l/J 

“« = CihpelJM'j - 1 ) 2 ) - UV]/J 

641 = v [(» 2 + An - 7P/W1 - 1) 2 )] 
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= iih p4/(pIi - !) 2 ) - v 2 \/ J 

and 

<12 = \/vev? 

Making these substitutions into equation (3.12) we obtain the reduced form 
of the algorithm 


/ + hS^UI + A c ) n + if) / + M„(V7 + B c ) n + Df 


A Q”- 
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Finally, to avoid the source term as well as to take advantage of existing 
codes, one can rewrite these equations in the form 


C 


I + hS({UI + A c r + Df> I + A«,(V7 + B c ) n + D& 


caq" 


— At 


S'E" + 6 V F* + D (4, <3 


(3.16) 

The right-hand-side of this equation is the same as equation (3.3a), while 

A 

the left-hand-side maintains the update for Q quantities rather than Q. The 


advantage of this last form is that one can readily take an existing code 
for equation (3.2) and modify only the left-hand-side operators so as to 
take advantage of the reduced inversion work. The steady state solution of 
equation (3.16) is clearly identical to that of equation (3.3a). 
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IV. Results 

The algorithm given by equation (3.16) was tested on a NACA 0012 airfoil 
under transonic flow conditions in order to verify that no adverse stability 
effects are incurred by using the matrix reductions. Versions of the basic 
NASA Ames Research Center AIR2D/ARC2D, which solve equation (3.3a), 
were modified to the form of equation (3.16). As noted earlier, only that 
portion of the code which deals with the left-hand-side operator has to be 
altered in order to implement the new algorithm. The basic code is described 
in refs. 5 and 6. 

A “C-type” topology was used for the airfoil calculations. The grids were 
generated with either a hyperbolic grid solver (ref. 12 and 13) that imposes 
orthogonality up to numerical errors of truncation and smoothing, or an 
algebraic construction that uses a method similar to the control function 
approach of Eiseman (ref. 14). In all cases, the far-field boundaries were 
placed approximately 20 chords from the body. Boundary conditions on the 
body, wake, and freestream are further described in refs. 5 and 6. 

As a first calculation, a relatively coarse grid was used which was generated 
with the hyperbolic grid solver. This grid has 157 points in the ^-direction 
and 33 points in the ^-direction (fig. 2a and 2b). The pressure distribution 
on the NACA 0012 for M ^ — 0.8 and a = 0. is shown in Figure 2c. 
This distribution was computed using both the standard algorithm, equation 
(3.3a), and the reduced matrix form, equation (3.16). Both algorithms were 
run using Euler implicit differencing, which is first-order accurate in time. 
As shown in fig. 2d, the residual histories are virtually identical. However, 
the standard algorithm required 120 CPU seconds of CRAY X-MP time per 
1000 iterations and the new algorithm required 54 CPU seconds per 1000 
iterations. 

Numerical experimentation, in which At was adjusted over a wide range, 
has confirmed that the new algorithm is as stable as the standard algorithm. 
As an example of the robustness of the formulation using equation (3.16), 
a much finer grid was used to compute the previous case. Figures 3a and 
3b show a view of the grid with 249 points in the ^-direction and 50 points 
in the ^-direction. The solution on this grid (Figures 3c-3d) was computed 
using a spacially varying time step scaled by 1 ^ y y-j where J is the metric 

Jacobian. For this calculation, the drag coefficient was converged to 5 digits 
in 600 iterations or 68 seconds of CRAY X-MP time. Figure 3e shows the 
residual history for this calculation. 

Recently, Beam and Bailey have reported in private communication (see 
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also ref. 6) that significantly faster convergence can be obtained by incor- 
porating the fourth-order smoothing terms implicitly into the Beam and 
Warming algorithm. For the standard algorithm, imposing fourth-order dis- 
sipation implicitly necessitates 4x4 block pentadiagonal inversions which are 
relatively expensive. The reduced matrix algorithm, however, requires only 
scalar and 2x2 block pentadiagonal inversions which are much less costly. 
In fact, a special 2x2 block pentadiagonal inversion algorithm, coded for 
this algorithm, requires only 107 arithmetic operations per entry. 

In testing the sound speed-velocity splitting algorithm with implicit fourth 
order dissipation, a more advanced form of numerical filter was implemented. 
This filter uses both second- and fourth-order differences such that the second- 
order difference is dominant near shock waves. The spectral radius of the 
flux Jacobian matrices is used as a scaling to the smoothing coefficients in 
an attempt to mimic the dissipative nature of second-order flux split upwind 
schemes (ref. 15). Full details are given in ref. 6. 

As a test case, an algebraic grid was used to compute the flow field around 
a NACA 0012 airfoil at — 0.8 and a = 1.25 degrees. The grid has 
193 points in the ^-direction and 33 points in the ^-direction. The grid and 
flow field solution are shown in fig. 4a through 4d. The residual history (fig. 
4e) shows the rapid convergence. For this particular case, the CPU time on 
the CRAY X-MP was 93 seconds per 1000 iterations. This compares with 65 
seconds per 1000 iterations using the tridiagonal version of the new algorithm 
with constant smoothing coefficients which had much slower convergence. It 
should be noted that the smoothing filter used for the pentadiagonal version 
required 16 seconds more CPU time than the same smoothing with a constant 
coefficient filter. 

Our computational results verified that the reduced algorithm was every bit 
as numerically stable as the standard algorithm. This experience is similar to 
what was observed previously by Steger on an algorithm version that reduced 
to a 3 x 3 block and a scalar. Here, with reduction to the 2x2 block, 
just over a factor of two was saved in overall computer time by using the 
reduced algorithm (3.16) in place of equation (3.3a). The precise savings that 
can be obtained is, of course, machine and coding dependent and a larger 
improvement should be obtained in three dimensions. 

Y. Conclusions 

By substituting similar reducible matrices for the Jacobian matrices of 
the flux vectors, a substantial reduction has been achieved in the arith- 
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metic operations needed in solving approximate factorization implicit al- 
gorithms such as Beam- Warming. Numerical calculations indicate that the 
new similarity split matrices can apparently be used without any loss of 
numerical stability, although time accuracy can be degraded unless additional 
steps are taken. The new method can be readily retrofit within existing codes, 
and can likely be extended to viscous flow calculations as well. 
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Appendix A 

Extension of Sound Speed and Velocity Splitting to the Navier-Stokes Equations 

The sound speed and velocity splitting concept can be applied to the 
Navier-Stokes equations because of the special form of the viscous Jacobians. 
The Navier-Stokes equations in strong conservation law form can be 
represented as 

®tQ ~h d x E T* dyF — d x E vx -f- d x E vy -j~ d y F vx d y F vy (A.l) 


where Q, E, and F have their previous definitions. The viscous flux terms 
are given by 

0 


with 


E vx = 


vx 


(4/3)//u* 

fiv x 

L F VX 4 

0 

t*v x 

~(2/3)fiu x 

E vx4 


E. 


vy 


F 

i 1 \ 


vy 


0 

(2/3)//% 

HUy 

E V yi 

' 0 ‘ 

»Uy 

(4/3)//% 

L E vy 4 


E vx 4 = (///2)<9 x [(4/3)u 2 + V 2 ] + VpdxC* 


(A.2a) 


(A.2a) 


E vyi = fiviiy — (2/3)//m% 

F vx 4 = fiuv x — (2/3 )flVU x 
E V y 4 = {(i/2)dy[u 2 + (4/3)y 2 ] + ///? d y c 2 


and 

P = Pr“ 1 (7 — l) - " 1 


Beam and Warming have developed a general class of integration techniques 
for equation (A.l) which is illustrated here for Euler implicit time differencing 

as 


AQ“ - At[S, A£" + S y AF n - i.AC - J ,Af 

= -A t\5,E" + t,F n - l,(K, + El v ) - 5,(f + F%)] ' ' 
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where A Q n ~ Q n + 1 — Q n and 6 is a midpoint operator. In this algorithm the 
cross derivative flux terms E vy and F vx are treated explicitly. As discussed by 
Beam and Warming, this circumvents the need to include an implicit viscous 
Jacobain for these cross terms, and it also maintains unconditional stablility 
for the scalar model equation. 

The spacial flux terms have been locally linearized in time. In locally 
linearizing these terms Beam and Warming expand the viscous terms in a 
Taylor series using both the function and its derivative, as for example 


A E n = E n+1 — E n 


'vx 


L dQ J 


(Q n+1 — Q") + 


a 1 (Q n+1 - Q n ) 

(A.4) 

Alternately, these terms can be expanded in terms of Q alone, for example 


dE v . 

L dQx J 


AE VX 


dE , 


vx 


[ dQ 


(Q 


n -\- 1 


Q n ) 


(A.5) 


where is now a differential operator. This latter form will be used here. 
Typically p, is not linearized in time and the Prandtl number is treated as if 
it were a constant value as far as time linearization is concerned. For a term 
such as j.idyU, fi is frozen in time and Hod y u expanded such as 


p 0 d y {pu/p) = iiody[{poUo/po) — (Po«oM> 2 )(/> ~ P o) + Po l {pu — p 0 u 0 )} 
or in general 

d x E vx — d x \E vx \o -f- [iod x {A. vx \o(Q Qo))]q (A.6a) 

dyF vy = 9 y [Fvy|o + Hody{B vy \o(Q — <5o))] 0 (A.6b) 

with 
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0 

0 

0 - 


0 0 

0 

0 - 

— ( 4 / 3 )(u /p) 

(4/3 IP” 1 

0 

o 

, B v y 

— {u/p) p 4 
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- ^il i>x 2 
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(A.7) 


and 


«,i = -(4/3)[a 2 //>] - W-h\ + + (« 2 + » 2 )//>] 
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a *2 = [(4/3) - P\{ufp) 

0*3 = (1 — PKv/p) 

b yl = ~[u 2 /p] - (4/3)[w 2 /p] + P[—(e/p) + (u 2 + v 2 )/p] 
by 2 = (1 ~ /?)(«/p) 
fc s/3 = [(4/3)- p](v/p) 

With the introduction of the sound speed and velocity matrix substitutions 
for the convection Jacobian matrices, we can obt ain operators L x and L y that 
are as reducible as before because of the form of A vx and B vy . To illustrate 
this we will drop the numerical dissipation terras as before. Then equation 
(A.3) becomes 

LtLyLQ" = - At\6 x E " + 6,F" - + K„) - SPK* + K,)] 

(A.8) 

where 
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Here I is the identity matrix and a c , etc., are elements of A c etc. Keep in 
mind that the viscous terms have been second-order accurate differenced to 
maintain tridiagonal structure using midpoint operators 6. 

We can now examine the structure of L x for example, and show that the 
equations uncouple to two scalar and one 2x2 block tridiagonal. The first 
row of L x is clearly uncoupled and can be solved for as a scalar tridiagonal. 
Its solution is then used to update the right-hand-side (RHS) of rows 2, 3, and 
4. With the first column of the matrices brought to the RHS, the third row 
is seen to be uncoupled and can thus be solved as another scalar tridiagonal. 
The RHS terms of the second and fourth rows can now be updated and what 
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remains is a 2 x 2 block tridiagonal to be solved. The L y operator inverts in 
a similar fashion. 

Thus, the addition of the viscous terms in Cartesian coordinates does 
not prevent the decoupling technique that was successfully used for the 
inviscid gasdynamic equations. Preliminary work shows that, using the local 
transformation of Section III, the same decoupling may be possible for the 
Navier-Stokes equations in generalized coordinates. This matter is being 
further investigated. 
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Appendix B 

Stability of Sound Speed and Velocity Splitting 

In order to study the stability of the algorithm proposed in Section H, 
the one-dimensional Euler equations will be tested for linear stability using a 
Fourier (or matrix) analysis. 

The one-dimensional Euler equations in Cartesian coordinates are given by 

d t Q + d x Er= 0 (B.l) 


where 


Q 


' p' 


• pu - 

pu 

, E = 

pu 2 -f- p 

. e . 


u{e + p\ 


A frozen coefficient form of equation (B.l) is given by 


(B.2) 


d t Q+Ad x E = 0 


(B-3) 


where A is the true Jacobian, [dE/dQ], and * denotes that A is evaluated 
using a frozen value of Q, Q*. Implicit first-order- accurate differencing of 
equation (B.3) is given in delta form as 


/ + hA 6 X 


(<3 n+1 — Q") = —hA'S x Q n 


(B.4) 


In our sound speed-velocity splitting scheme, we replace the A* matrix of the 
left-hand-side operator with the similarity matrices u / -f- A c to obtain 


I-f- h{u! + A? c )5: 


(Q n+1 - Q n ) = -hjChQ" 


(B.5) 


The right-hand-side matrix A * is not altered in our procedure, and must not 
be changed in equation (B.5) as it represents the correct linearization of E 
about Q*. The question to be addressed is, is equation (B.5) as stable as 
equation (B.4) now that the implicit operator has been altered? 

If we were to perform a stability analysis of equation (B.4), we would 
diagonalize A* using its eigenvectors. We can then readily perform the 
stability analysis for difference equations corresponding to scalar partial 
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differential equations. For equation (B.5) we can use the eigenvectors of A * 
to diagonalize u I -J- A c , but the right-hand-side matrix will not be brought 
into diagonal form. Nevertheless, we find that the resulting equations can 
still be analyzed. 

The matrix A c for one-dimensional flow is given by 


A c = { 7 - 1) 


r o 

u 2 / 2 


o On 

— u 1 


ah °42 


u. 


(B.6) 


where 


a c u = u[u 2 / 2 - c 2 /( 7 - l) 2 ] 


a 42 = -u 2 + c 2 /( 7- l) 2 
Eigenvectors matrices of A c are given by 

- [(7 - 1 )u 2 + 2cu]/(4c) [(7 - 1 )u + c]/(2c) 


X" 


(7- l)/(2c)' 


[(7- l)u 2 - 2c«]/(4c) -[(7- l)u~ c]/(2c) (7— ])/(2c) 
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0 


0 
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1 1 


0 0 

1 1 u 

U - (c/(7 - 1)) u + (c/(7 - 1)) u 2 /2j 


(B.7) 

(B.8) 


and 




— e 0 O' 

0 c 0 (B.9) 

1 0 0 0. 

Multiplying equation (B.5) by the frozen-coefficient inverse of the eigenvec- 
tors (X*)" 1 gives 


I h{ U I ~b 


(x , r 1 (Q" +1 ~Q n ) = -h{x*r'A*x%(x*r 1 Q n 


or in nondelta form with W — (X*) 


(B.10) 


I -j- h(u I -j- 


W n+1 ^ I-h[(X*)^ 1 A*X*-uI-A\]6 x W n (B.ll) 
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However, (AT*) *A*X* — u* I — 


is a very simple matrix 


(X*)~ 1 A*X* - u*I - A* Ac 


•0 0 O' 

0 0 0 

.1 1 0. 


(B.12) 


Assuming periodic boundary conditions so as to use a Fourier stability 
analysis (or the matrix stability method using circulant matrices), the finite 
difference operator S x can be transformed to 

2 isin(9j) 

- — = K 

2Ax 


and W denotes Fourier transformation of W ( i.e. transformation via the 
eigenvectors of a circulant matrix). 

Applying the Fourier method to equation (B. 11) and taking the inverse of 
the left hand side operator we obtain 


(B. 13) 


(B.14) 


The amplification matrix associated with equation (B.14) is lower triangular. 
Therefore, its eigenvalues are the diagonal elements, and these clearly have 
modulus less than one for any positive value of h — At. In fact, these 
eigenvalues are identical to those obtained with the standard scheme given 
by equation (B.4). Thus the necessary condition for stability is met for any 
values of h. A sufficient condition for stability requires bounding the norm of 
the 3x3 amplification matrix. If we observe the 3, 1 and 3, 2 elements, it is 
clear that the £<& norm can be unbounded for u -+ 0 and h -*• oo; however, 
even in this norm, powers of the amplification matrix will be quickly bounded 
unless u is identically zero. 

It is interesting to note that if one were to form A c using the Jacobian of 
E Ci the resulting algorithm is unconditionally unstable. In this case A c is 


or 
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given as 


where 


A c = (q ■— 1 ) 


r 0 

u 2 / 2 

L a 41 
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■u 
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42 


0- 
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. I 

u_ 


(B.15) 


a\ i = u\u 2 /2- c 2 / 7 ( 7 - 1 )] 

« 42--^ 2 + c 2 / 7 ( 7 - 1 ) 

If we procede in a fashion similar to the above, the resultant representation 
in Fourier space for equation (B.5) is 
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— ihK 
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a — \Ay(7 

-i) /? = v^ 




The maximum eigenvalue of the amplification matrix in equation (B.16) 
can be shown to be greater than one for all At’s greater than zero, and 
the necessary condition for stability cannot be met. In actual numerical 
experimentation where numerical dissipation is added, small values of At 
could be used to obtain stable results. 
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SUBBLOCK SIZE, M 

Figure 1. - Arithmetic operation counts of block-tridiagonals using L-U decomposition. 
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(a) (b) 



Figure 2. - Coarse grid solution about NACA 0012 airfoil using reduced 
matrix algorithm, (a) Overview of 157 X 33 grid generated using hyperbolic 
grid generator, (b) Detail near body, (c) Pressure distribution = 0.8 
and a = 0.0). (d) Residual history comparison between reduced matrix and 
standard algorithms. 


28 







Figure 3. - Fine grid solution about NACA 0012 airfoil using reduced matrix 
algorithm with scaled At. (a) Overview of 249 X 50 grid generated using 
hyperbolic grid generator, (b) Detail near body, (c) Pressure distribution on 
airfoil (M^ = 0.8 and a — 0.0). (d) Pressure contours near airfoil. 
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(e) Residual history for reduced matrix algorithm with scaled At. 

Figure 3. - Concluded. 
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Figure 4. - Solution about NACA 0012 airfoil using advanced numerical 
filter, scaled A t, and pentadiagonal inversions, (a) Overview of 193 x 33 grid 
generated using algebraic grid generator, (b) Detail near body.(c) Pressure 
distribution on airfoil (M^ = 0.8 and a = 1.25). (d) Residual history for 
reduced matrix algorithm. 








(e) Residual history for reduced matrix algorithm with scaled 
At and pentadiagonal inversions. 

Figure 4. - Concluded. 
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